This workflow is based on using the R package kwb.hydrus1d, which is compatible with the last official release of HYDRUS1D v4.17.0140. The development of this new R package was necessary, due to the fact that already available model wrappers were developed for outdated HYDRUS1D versions (e.g. python package phydrus requires HYDRUS1D v4.08) and/or are not well maintained (e.g. open issues/pull requests for R package hydrusR).
The workflow below describes all the steps required for:
Input data preparation
Modifying (i.e.
atmosphericinput data defined in fileATMOSPHERE.in) an already pre-preparedHYDRUS1Dmodel (using theHYDRUS1D GUI)Running it from within
Rby using the command line (cmd) andImporting and analysing the
HYDRUS1Dresults within R.
Install R Package
# Enable this universe
options(repos = c(
kwbr = 'https://kwb-r.r-universe.dev',
CRAN = 'https://cloud.r-project.org'))
# Install R package
install.packages('flextreat.hydrus1d')Define Paths
paths_list <- list(
extdata = system.file("extdata", package = "flextreat.hydrus1d"),
exe_dir = "<extdata>/model",
model_name = "test",
model_dir = "<exe_dir>/<model_name>",
atmosphere = "<model_dir>/ATMOSPH.IN",
a_level = "<model_dir>/A_LEVEL.out",
t_level = "<model_dir>/T_LEVEL.out",
runinf = "<model_dir>/Run_Inf.out",
solute_id = 1,
solute = "<model_dir>/solute<solute_id>.out",
soil_data = "<extdata>/input-data/soil/soil_geolog.csv"
)
paths <- kwb.utils::resolve(paths_list)Input Data
Soil
Soil data is derived from depth-dependent grain-size analysis of soil samples taken in Braunschweig. The following required input parameters for the van Genuchten model used in HYDRUS1D were derived by using the pedotransfer function ROSETTA-API version 3, which was developed by Zhang and Schaap, 2017.
library(flextreat.hydrus1d)
soil_dat <- readr::read_csv(paths$soil_data, show_col_types = FALSE) %>%
dplyr::mutate(layer_id = dplyr::if_else(.data$tiefe_cm_ende < 60,
1,
2),
thickness_cm = .data$tiefe_cm_ende - .data$tiefe_cm_start,
sand = .data$S_prozent + .data$G_prozent,
silt = .data$U_prozent,
clay = round(100 - .data$sand - .data$silt,
1)) %>%
dplyr::mutate(clay = dplyr::if_else(.data$clay < 0,
0,
.data$clay)) %>%
soilDB:::ROSETTA(vars = c("sand", "silt", "clay"),
v = "3") %>%
dplyr::mutate(alpha = 10^.data$alpha,
npar = 10^.data$npar,
ksat = 10^.data$ksat)
knitr::kable(soil_dat)| tiefe_cm | tiefe_cm_start | tiefe_cm_ende | gluehverlust_prozent | U_prozent | S_prozent | G_prozent | kf_beyer | layer_id | thickness_cm | sand | silt | clay | theta_r | theta_s | alpha | npar | ksat | .rosetta.model | .rosetta.version |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 - 30 | 0 | 30 | 2.7 | 6.4 | 93.1 | 0.4 | 9.5e-05 | 1 | 30 | 93.5 | 6.4 | 0.1 | 0.0481310 | 0.3671023 | 0.0334257 | 3.014652 | 716.3720 | 2 | 3 |
| 30 - 55 | 30 | 55 | 1.3 | 6.3 | 93.3 | 0.4 | 1.3e-04 | 1 | 25 | 93.7 | 6.3 | 0.0 | 0.0480175 | 0.3669818 | 0.0335768 | 3.053357 | 739.3399 | 2 | 3 |
| 55 - 80 | 55 | 80 | 0.6 | 1.7 | 97.3 | 1.1 | 1.9e-04 | 2 | 25 | 98.4 | 1.7 | 0.0 | 0.0485529 | 0.3626158 | 0.0355315 | 3.951321 | 1296.6154 | 2 | 3 |
| 80 - 150 | 80 | 150 | 0.2 | 1.3 | 98.6 | 0.1 | 2.0e-04 | 2 | 70 | 98.7 | 1.3 | 0.0 | 0.0486320 | 0.3621635 | 0.0356286 | 4.021821 | 1343.1455 | 2 | 3 |
| 150 - 190 | 150 | 190 | 0.2 | 1.0 | 98.7 | 0.3 | 2.2e-04 | 2 | 40 | 99.0 | 1.0 | 0.0 | 0.0486788 | 0.3618558 | 0.0357316 | 4.088764 | 1388.8655 | 2 | 3 |
| 190 - 210 | 190 | 210 | 0.3 | 2.5 | 97.2 | 0.3 | 1.7e-04 | 2 | 20 | 97.5 | 2.5 | 0.0 | 0.0484589 | 0.3633709 | 0.0351930 | 3.764078 | 1171.5976 | 2 | 3 |
Due to similar soil characteristics, only two different layers (column layer_id) were defined: - layer 1: ranging from 0 cm - 55 cm depth - layer 2: ranging from 55 cm - 210 cm depth
The spatially aggregation for each layer (geometric mean) of the input parameters for the van Genuchten model is performed with the code below and shown in the subsequent table.
soil_dat_per_layer <- soil_dat %>%
dplyr::group_by(.data$layer_id) %>%
dplyr::summarise(layer_thickness = max(.data$tiefe_cm_ende) - min(.data$tiefe_cm_start),
theta_r = sum(.data$theta_r*.data$thickness_cm)/.data$layer_thickness,
theta_s = sum(.data$theta_s*.data$thickness_cm)/.data$layer_thickness,
alpha = sum(.data$alpha*.data$thickness_cm)/.data$layer_thickness,
npar = sum(.data$npar*.data$thickness_cm)/.data$layer_thickness,
ksat = sum(.data$ksat*.data$thickness_cm)/.data$layer_thickness)
knitr::kable(soil_dat_per_layer)| layer_id | layer_thickness | theta_r | theta_s | alpha | npar | ksat |
|---|---|---|---|---|---|---|
| 1 | 55 | 0.0480794 | 0.3670475 | 0.0334944 | 3.032245 | 726.812 |
| 2 | 155 | 0.0486090 | 0.3623128 | 0.0355833 | 3.994469 | 1325.304 |
Atmospheric Boundary Conditions
In total three different atmospheric input boundary conditions (precipitation, potential evaporation and irrigation) are used within the HYDRUS1D model and described in the following subchapters in more detail.
Rainfall
Rainfall is based on historical hourly raw data downloaded from DWD open-data ftp server for station_id = 662 (i.e. Braunschweig). Data is aggregated within R to daily sums by using the code below:
install.packages("rdwd")
library(dplyr)
rdwd::updateRdwd()
rdwd::findID("Braunschweig")
rdwd::selectDWD(name = "Braunschweig", res = "daily")
url_bs_rain <- rdwd::selectDWD(name = "Braunschweig",
res = "hourly",
var = "precipitation",
per = "historical" )
bs_rain <- rdwd::dataDWD(url_bs_rain)
precipitation_hourly <- rdwd::dataDWD(url_bs_rain) %>%
dplyr::select(.data$MESS_DATUM, .data$R1) %>%
dplyr::rename("datetime" = "MESS_DATUM",
"precipitation_mm" = "R1")
usethis::use_data(precipitation_hourly)
precipitation_daily <- precipitation_hourly %>%
dplyr::mutate("date" = as.Date(datetime)) %>%
dplyr::group_by(date) %>%
dplyr::summarise(rain_mm = sum(precipitation_mm))
usethis::use_data(precipitation_daily)
DT::datatable(head(flextreat.hydrus1d::precipitation_daily))Potential Evaporation
Is based on daily potential evaporation grids (1km x 1km) downloaded from DWD open-data server, where all grids (i.e. 46) within the Abwasserverregnungsgebiet.shp were selected and spatially aggregated (mean, sd, min, max) by using the code below:
shape_file <- system.file("extdata/input-data/gis/Abwasserverregnungsgebiet.shp",
package = "flextreat.hydrus1d")
# Only data of full months can currently be read!
evapo_p <- kwb.dwd::read_daily_data_over_shape(
file = shape_file,
variable = "evapo_p",
from = "201701",
to = "202012"
)
usethis::use_data(evapo_p)Irrigation
Monthly irrigation volumes were provided by Abwasserverband Braunschweig for the time period 2017-01-01 - 2020-12-31 as csv file and separates between two sources (groundwater and clearwater). Data preparation was carried out with the code below:
irrigation_file <- system.file("extdata/input-data/Beregnungsmengen_AVB.csv",
package = "flextreat.hydrus1d")
irrigation_area <- rgdal::readOGR(dsn = shape_file)
irrigation <- read.csv2(irrigation_file) %>%
dplyr::select(- .data$Monat) %>%
dplyr::rename(irrigation_m3 = .data$Menge_m3,
source = .data$Typ,
month = .data$Monat_num,
year = .data$Jahr) %>%
dplyr::mutate(date_start = as.Date(sprintf("%d-%02d-01",
.data$year,
.data$month)),
days_in_month = as.numeric(lubridate::days_in_month(.data$date_start)),
date_end = as.Date(sprintf("%d-%02d-%02d",
.data$year,
.data$month,
.data$days_in_month)),
source = kwb.utils::multiSubstitute(.data$source,
replacements = list("Grundwasser" = "groundwater.mmPerDay",
"Klarwasser" = "clearwater.mmPerDay")),
irrigation_cbmPerDay = .data$irrigation_m3/.data$days_in_month,
irrigation_area_sqm = irrigation_area$area,
irrigation_mmPerDay = 1000*irrigation_cbmPerDay/irrigation_area_sqm) %>%
dplyr::select(.data$year,
.data$month,
.data$days_in_month,
.data$date_start,
.data$date_end,
.data$source,
.data$irrigation_mmPerDay,
.data$irrigation_area_sqm) %>%
tidyr::pivot_wider(names_from = .data$source,
values_from = .data$irrigation_mmPerDay)
usethis::use_data(irrigation)
DT::datatable(head(flextreat.hydrus1d::irrigation))HYDRUS-1D
Modify Input Files
All model input files were initially setup by using the HYDRUS1D-GUI but the following below were modified manually with the
SELECTOR.in
Soil input data were entered manually in SELECTOR.in for the two layers defined the second table in Chapter: Input Data - Soil.
ATMOSPHERE.in
Based on the atmospheric input data (see Chapter: Atmospheric Boundary Conditions) the ATMOSPHERE.in file for HYDRUS1D is prepared with the code below and starts with the hydrological summer half year (assumption: soil is fully wetted at end of winter half year):
atm <- flextreat.hydrus1d::prepare_atmosphere_data()
atm_selected <- flextreat.hydrus1d::select_hydrologic_years(atm)
atm_prep <- flextreat.hydrus1d::prepare_atmosphere(
atm_selected,
conc_irrig_clearwater = 1, # use as "clearwater" tracer
conc_irrig_groundwater = 0,
conc_rain = 0)
atmos <- kwb.hydrus1d::write_atmosphere(
atm = atm_prep,
MaxAL = nrow(atm_prep),
round_digits = 6 # increase precision for "clearwater" tracer!
)
writeLines(atmos, paths$atmosphere)The selected time period covers 1278 days (2017-05-01 - 2020-10-30), i.e. it covers 7 hydrological half years.
DT::datatable(atm_selected)These time-series were converted with the function flextreat.hydrus1d::prepare_atmosphere() into the data format required by HYDRUS1D. Due to the fact, that irrigation rates (i.e. sum of clearwater.mmPerDay and groundwater.mmPerDay) cannot be entered separately as input column within HYDRUS1D, these were simply added to th prec (i.e. precipitation) column. The whole time series defined in ATMOSPHERE.in is shown below:
DT::datatable(atm_prep)Run Model
Finally the model is run automatically by using the following code:
exe_path <- kwb.hydrus1d::check_hydrus_exe(dir = paths$exe_dir,
skip_preinstalled = TRUE)
#> Checking if download of HYDRUS1D executable v4.17.0140 from 'https://github.com/mrustl/hydrus1d/archive/refs/tags/v4.17.0140.zip' was successful ... ok. (0.00s)
kwb.hydrus1d:::run_model(exe_path = exe_path,
model_path = paths$model_dir)
#> Warning in shell(cmd = sprintf("cd %s && %s", fs::path_abs(target_dir), : 'cd
#> D:/a/_temp/Library/flextreat.hydrus1d/extdata/model && H1D_CALC.exe' execution
#> failed with error code 24Read Results
runinf <- kwb.hydrus1d::read_runinf(paths$runinf)
summary(runinf)
#> t_level time dt itr_w
#> Min. : 1 Min. : 0.05 Min. :0.0001065 Min. : 2.000
#> 1st Qu.: 5954 1st Qu.: 285.35 1st Qu.:0.0220900 1st Qu.: 2.000
#> Median :11908 Median : 649.70 Median :0.0495000 Median : 2.000
#> Mean :11908 Mean : 624.94 Mean :0.0537056 Mean : 2.483
#> 3rd Qu.:17862 3rd Qu.: 957.35 3rd Qu.:0.0944970 3rd Qu.: 2.000
#> Max. :23815 Max. :1279.00 Max. :0.1000000 Max. :20.000
#> itr_c it_cum kod_t kod_b converg
#> Min. :1 Min. : 4 Min. :-4.000000 Min. :-5 Mode:logical
#> 1st Qu.:1 1st Qu.:19307 1st Qu.:-4.000000 1st Qu.:-5 TRUE:23815
#> Median :1 Median :36910 Median : 4.000000 Median :-5
#> Mean :1 Mean :37823 Mean : 0.004199 Mean :-5
#> 3rd Qu.:1 3rd Qu.:56331 3rd Qu.: 4.000000 3rd Qu.:-5
#> Max. :1 Max. :75157 Max. : 4.000000 Max. :-5
#> peclet courant
#> Min. :0.1 Min. :0.0000
#> 1st Qu.:0.1 1st Qu.:0.0150
#> Median :0.1 Median :0.0400
#> Mean :0.1 Mean :0.0829
#> 3rd Qu.:0.1 3rd Qu.:0.1170
#> Max. :0.1 Max. :0.9960
a_level <- kwb.hydrus1d::read_alevel(paths$a_level)
a_level
#> # A tibble: 1,279 × 10
#> time sum_r_top sum_r…¹ sum_v…² sum_v…³ sum_v_…⁴ h_top h_root h_bot a_level
#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1 0.246 0 0.0672 0 -0.00173 -1 e5 0 -116. 1
#> 2 2 -0.213 0 -0.392 0 -0.00346 -8.30e1 0 -116. 2
#> 3 3 -0.0375 0 -0.216 0 -0.00520 -1.67e3 0 -116. 3
#> 4 4 -0.162 0 -0.341 0 -0.00693 -1.03e2 0 -116. 4
#> 5 5 -0.157 0 -0.335 0 -0.00866 -1.23e2 0 -116. 5
#> 6 6 -0.0311 0 -0.221 0 -0.0104 -1 e5 0 -116. 6
#> 7 7 0.188 0 -0.187 0 -0.0121 -1 e5 0 -116. 7
#> 8 8 0.232 0 -0.167 0 -0.0139 -1 e5 0 -116. 8
#> 9 9 0.408 0 -0.151 0 -0.0156 -1 e5 0 -116. 9
#> 10 10 0.583 0 -0.139 0 -0.0173 -1 e5 0 -116. 10
#> # … with 1,269 more rows, and abbreviated variable names ¹sum_r_root,
#> # ²sum_v_top, ³sum_v_root, ⁴sum_v_bot
t_level <- kwb.hydrus1d::read_tlevel(paths$t_level)
t_level
#> # A tibble: 23,815 × 22
#> time r_top r_root v_top v_root v_bot sum_r_top sum_r_r…¹ sum_v…² sum_v…³
#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.05 0.246 0 0.246 0 -0.00173 0.0123 0 0.0123 0
#> 2 0.1 0.246 0 0.246 0 -0.00173 0.0246 0 0.0246 0
#> 3 0.15 0.246 0 0.246 0 -0.00173 0.0369 0 0.0369 0
#> 4 0.168 0.246 0 0.246 0 -0.00173 0.0414 0 0.0414 0
#> 5 0.188 0.246 0 0.0664 0 -0.00173 0.0464 0 0.0427 0
#> 6 0.203 0.246 0 0.0433 0 -0.00173 0.0498 0 0.0433 0
#> 7 0.218 0.246 0 0.0562 0 -0.00173 0.0536 0 0.0442 0
#> 8 0.235 0.246 0 0.0596 0 -0.00173 0.0578 0 0.0452 0
#> 9 0.254 0.246 0 0.0575 0 -0.00173 0.0625 0 0.0463 0
#> 10 0.275 0.246 0 0.0494 0 -0.00173 0.0675 0 0.0473 0
#> # … with 23,805 more rows, 12 more variables: sum_v_bot <dbl>, h_top <dbl>,
#> # h_root <dbl>, h_bot <dbl>, run_off <dbl>, sum_run_off <dbl>, volume <dbl>,
#> # sum_infil <dbl>, sum_evap <dbl>, t_level <dbl>, cum_w_trans <dbl>,
#> # snow_layer <dbl>, and abbreviated variable names ¹sum_r_root, ²sum_v_top,
#> # ³sum_v_root
solute <- kwb.hydrus1d::read_solute(paths$solute)
solute
#> # A tibble: 23,815 × 19
#> time cv_top cv_bot sum_cv…¹ sum_c…² cv_ch0 cv_ch1 c_top c_root c_bot cv_root
#> * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 0.05 1.57 0 0.0783 0 0 0 1.00 0 0 0
#> 2 0.1 0.779 0 0.117 0 0 0 1.00 0 0 0
#> 3 0.15 0.558 0 0.145 0 0 0 1.00 0 0 0
#> 4 0.168 0.444 0 0.153 0 0 0 1.00 0 0 0
#> 5 0.188 0.298 0 0.159 0 0 0 1.00 0 0 0
#> 6 0.203 0.143 0 0.161 0 0 0 1.00 0 0 0
#> 7 0.218 0.129 0 0.163 0 0 0 1.00 0 0 0
#> 8 0.235 0.138 0 0.166 0 0 0 1.00 0 0 0
#> 9 0.254 0.134 0 0.168 0 0 0 1.00 0 0 0
#> 10 0.275 0.120 0 0.171 0 0 0 1.00 0 0 0
#> # … with 23,805 more rows, 8 more variables: sum_cv_root <dbl>,
#> # sum_cv_n_eql <dbl>, t_level <dbl>, c_gwl <dbl>, c_run_off <dbl>,
#> # sum_c_run_off <dbl>, cv_i <dbl>, sum_cv_i <dbl>, and abbreviated variable
#> # names ¹sum_cv_top, ²sum_cv_botPlot
gg <- solute %>%
dplyr::select(.data$time, .data$c_top, .data$c_bot) %>%
dplyr::mutate(date = as.POSIXct("2017-05-01", tz = "UTC") + .data$time*24*3600) %>%
tidyr::pivot_longer(cols = c("c_top", "c_bot"),
names_to = "parameter",
values_to = "value") %>%
ggplot2::ggplot(mapping = ggplot2::aes(x = .data$date,
y = .data$value,
col = .data$parameter )) +
ggplot2::geom_point() +
ggplot2::geom_line() +
ggplot2::ylim(c(0,1)) +
ggplot2::theme_bw()
gg
plotly::ggplotly(gg)